function [P_plan,P_real,V,delta_T,fatigue,T0] = get_P(GR,CP,W,sigma,hight,V_wind,theta,Humidity,Temperature,turn_sign)
mT = 85;%总质量
fatigue = zeros(1,length(GR)+1); %疲劳因子
tensity = 1; %GR影响规划功率的系数因子
P_real = ones(1,length(GR));%预测的实际功率
P_plan = ones(1,length(GR));%规划的功率
V = ones(1,length(GR));%速度
v=0;%急弯上的切向速度
u=zeros(1,length(GR));%急弯上的向心速度
decrease =0.925;%入急弯速度衰减
r=4;%急弯半径
delta_T = zeros(1,length(GR));
g = 9.80665; %the acceleration of gravity at sea level
C_RR = 0.0032;%滚动阻力系数（包括轮胎和地面特性）

beta=(C_RR*mT*g);
p0 = 101325;
L=0.0065;
T0=288.15;
M = 0.0289644;
R_g = 8.31447;
R_d = 287.058;
R_v = 461.495;


for i=1:1:length(GR)

    rho = p0.*((1-L.*hight(1,i)./T0).^(g.*M./(R_g.*L))) ./(R_d.*Temperature) +...
    0.061078.* Humidity.*(1./(R_v.*Temperature) - 1./(R_d.*Temperature)).*...
    (10.^((7.5.*Temperature-2048.625)./(Temperature-35.85)));     
    
    alpha=1/2*rho*(0.2565+0.0044);

    if(abs(GR(1,i))<=0.01)
       P_plan(1,i) = CP;
    end
    if(GR(1,i)>0.01 && GR(1,i)<=0.03)
       P_plan(1,i) = CP.*(1+0.7.*tensity.*GR(1,i));
    end
    if(GR(1,i)>0.03 && GR(1,i)<=0.05)
       P_plan(1,i) = CP.*(1+0.7.*tensity.*GR(1,i));
    end
    if(GR(1,i)>0.05 && GR(1,i)<=0.1)
       P_plan(1,i) = CP.*(1 + tensity.*0.3) + ((GR(1,i)-0.04).^2) + 0.5*sigma*sin(10.*GR(1,i));
    end
    if(GR(1,i)<-0.01)
       P_plan(1,i) = CP.*0.7 + exp((GR(1,i)+0.01));
    end
    if(GR(1,i)>0.1)
       P_plan(1,i) = CP.*(1 + tensity);
    end

    
    P_real(1,i) = P_plan(1,i)*(1-fatigue(1,i));

    temp=Solve3Polynomial(alpha,V_wind(1,i).*alpha.*2,(alpha.*(V_wind(1,i).^2)+beta*cos(atan(GR(1,i)))+g.*mT*sin(atan(GR(1,i)))),-P_real(1,i));

    if(turn_sign(1,i)==0)
       V(1,i)=max(temp(imag(temp)==0));
    end

    if(turn_sign(1,i)==1)
       v =max(temp(imag(temp)==0)).*decrease;
       u(1,i)=4.*(100.*pi./12./85./abs((v-V(1,i-1))./delta_T(1,i-1))-1);
       V(1,i)=sqrt(v.^2 + u(1,i).^2);
    end

    fatigue(1,i+1) = fatigue(1,i) + 0.00001.*atan(P_real(1,i) - CP);
    
    w_real = 0;
    w_real = w_real + (P_real(1,i) - CP)* delta_T(1,i);
    if(w_real>W)
        fatigue(1,i+1) = 3.*(fatigue(1,i) + 0.00001.*atan(P_real(1,i) - CP));
    end

    if fatigue(1,i+1) <= 0
       fatigue(1,i+1) = 0;
    end
    

    if(i~=length(GR))
    delta_T(1,i)=1./((V(1,i)+V(1,i+1))./2);
    end
    if(i==length(GR))
    delta_T(1,i)=1./V(1,length(GR));
    end

end

t(1,1) = 1./((V(1,1)+V(1,2))./2);
for i=1:1:length(GR)-1
    t(1,i+1)=t(1,i)+1./((V(1,i)+V(1,i+1))./2);
end
t(1,length(GR))=1./V(1,length(GR))+t(1,length(GR)-1);
T0 = t(1,length(GR))./60;

%      T0 = sum(delta_T)/60;

if(w_real>CP.*T0)
    error("w_real is bigger than CP")
end


end

